
#INEQUALITY PLOTS
individual_level <-
  individual_level %>%
  mutate(ineq_level = case_when(ineq < median(ineq, na.rm = T)  ~ "low",
                                ineq > median(ineq, na.rm = T)  ~ "high",
                                TRUE ~ NA_character_))
low_ineq =
  individual_level %>%
  filter(dist < 10000) %>%
  filter(ineq_level == "low")

high_ineq =
  individual_level %>%
  filter(dist < 10000) %>%
  filter(ineq_level == "high")

model1_low_ineq <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing, data=low_ineq)
model2_low_ineq <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing, data=low_ineq)
model3_low_ineq <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing, data=low_ineq)
model4_low_ineq <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing, data=low_ineq)

model1_high_ineq <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing, data=high_ineq)
model2_high_ineq <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing, data=high_ineq)
model3_high_ineq <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing, data=high_ineq)
model4_high_ineq <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing, data=high_ineq)

model1_low_ineq <- coeftest(model1_low_ineq, vcov=cluster.vcov(model1_low_ineq, cluster = low_ineq$mkid14))
model2_low_ineq <- coeftest(model2_low_ineq, vcov=cluster.vcov(model2_low_ineq, cluster = low_ineq$mkid14))
model3_low_ineq <- coeftest(model3_low_ineq, vcov=cluster.vcov(model3_low_ineq, cluster = low_ineq$mkid14))
model4_low_ineq <- coeftest(model4_low_ineq, vcov=cluster.vcov(model4_low_ineq, cluster = low_ineq$mkid14))

model1_high_ineq <- coeftest(model1_high_ineq, vcov=cluster.vcov(model1_high_ineq, cluster = high_ineq$mkid14))
model2_high_ineq <- coeftest(model2_high_ineq, vcov=cluster.vcov(model2_high_ineq, cluster = high_ineq$mkid14))
model3_high_ineq <- coeftest(model3_high_ineq, vcov=cluster.vcov(model3_high_ineq, cluster = high_ineq$mkid14))
model4_high_ineq <- coeftest(model4_high_ineq, vcov=cluster.vcov(model4_high_ineq, cluster = high_ineq$mkid14))

model1_low_ineq <- tidy(model1_low_ineq) %>% mutate(outcome = "outcome1", level = "low")
model2_low_ineq <- tidy(model2_low_ineq) %>% mutate(outcome = "outcome2", level = "low")
model3_low_ineq <- tidy(model3_low_ineq) %>% mutate(outcome = "outcome3", level = "low")
model4_low_ineq <- tidy(model4_low_ineq) %>% mutate(outcome = "outcome4", level = "low")

model1_high_ineq <- tidy(model1_high_ineq) %>% mutate(outcome = "outcome1", level = "high")
model2_high_ineq <- tidy(model2_high_ineq) %>% mutate(outcome = "outcome2", level = "high")
model3_high_ineq <- tidy(model3_high_ineq) %>% mutate(outcome = "outcome3", level = "high")
model4_high_ineq <- tidy(model4_high_ineq) %>% mutate(outcome = "outcome4", level = "high")
  
inequality_plot <-
  bind_rows(model1_low_ineq, model2_low_ineq, model3_low_ineq, model4_low_ineq,
            model1_high_ineq, model2_high_ineq, model3_high_ineq, model4_high_ineq) %>%
  mutate(outcome_label = case_when(outcome == "outcome1" ~ "Place of worship",
                                   outcome == "outcome2" ~ "Live in village",
                                   outcome == "outcome3" ~ "Live in neighborhood",
                                   outcome == "outcome4" ~ "Rent room",
                                   TRUE ~ NA_character_)) %>%
  filter(term == "opium_legal") %>%
  mutate(Inequality = level) %>%
  ggplot(aes(x=estimate, y = outcome_label, group = Inequality, shape = Inequality, color = Inequality)) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbarh(aes(xmin=estimate-1.67*std.error, xmax = estimate+1.67*std.error), position = position_dodge(width = 0.4), height = 0) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  xlab("Effect of Opium Concession System") +
  theme_bw() +
  scale_color_grey() +
  theme(axis.ticks.y.left = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title.y = element_blank(),
        panel.grid.major.x = element_blank(),
        #legend.title = element_blank(),
        axis.line.y.left = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))


ggsave("./_4_outputs/figures/figure_8.tiff", plot = inequality_plot,  width = 5, height = 3, dpi = 300, compression = "lzw")
